Até o momento, vocês tiveram contato com uma série de pacotes e funções do R, incluindo pacotes próprios para organização de dataframes (dplyr) e visualização de dados (ggplot), além de funções para descrição e análise.
Hoje, vamos combinar e aprofundar as habilidades que vocês desenvolveram até aqui.
Manipular um grande conjunto de dados reais
Analisar um grande conjunto de dados reais
Discutir questões teóricas envolvendo múltiplas comparações de uma mesma VD
“[P]rojeto de tecnologia que usa inteligência artificial para auditar contas públicas e auxiliar no controle social”
Para a aula de hoje: usaremos o conjunto de dados camara, organizado pelo Prof Marcus Nunes (DEST/UFRN) a partir dos dados auditados pelo Serenata de Amor.
Os dados consistem de reembolsos pagos aos deputados a partir de sua quota parlamentar, que prevê gastos com refeições (apenas para os deputados), combustíveis, publicidade, passagens aéreas etc.
Vamos analisar um pequeno conjunto de dados, mas se você quiser baixar o conjunto todo, dê uma lida no blog do Prof. Marcus para saber como fazê-lo.
O conjunto reúne todos os reembolsos dados a todos os deputados entre 2009 e 2017.
No total são 1188 deputados, que receberam reembolsos por 22 tipos diferentes de subquotas.
Hoje, vamos analisar os dados em camara_sample.csv. Esse conjunto selecionou aleatoriamente um único gasto por deputados em cada uma das 22 subquotas.
Por que selecionar uma amostra? Porque o teste-t pressupõe independência dos dados, lembra?
Precisamos que haja apenas UM dado de cada deputado para cada tipo de subquota. Se houver mais de uma medida por deputado, precisaríamos usar outros teste (como modelos lineares mistos).
Importe o conjunto de dados camara_sample.csv usando read.csv.
Salve seu novo dataframe com o nome camara.sample.
camara.sample = read.csv("~/Documents/Mahayana/academics/cursos/introstatistics/Aulas/Aula7/camara_sample.csv")Nome do deputado (congressperson_name)
Estado representado (state)
Partido político (party)
Valor de reembolso (total_net_value)
Descrição da subquota (subquota_description)
Prestador de serviço (supplier)
Ano da despesa (year)
Seria bom conhecermos os níveis de cada uma das variáveis da nossa tabela. A função unique() faz isso. No nosso caso, sabermos os níveis de subquota_description será muito útil.
unique(camara.sample$subquota_description)
unique(camara.sample$state)unique(camara.sample$subquota_description)| x |
|---|
| Automotive vehicle renting or watercraft charter |
| Congressperson meal |
| Consultancy, research and technical work |
| Flight ticket issue |
| Flight tickets |
| Fuels and lubricants |
| Locomotion, meal and lodging |
| Lodging, except for congressperson from Distrito Federal |
| Maintenance of office supporting parliamentary activity |
| Postal services |
| Publication subscriptions |
| Publicity of parliamentary activity |
| Purchase of office supplies |
| Security service provided by specialized company |
| Software purchase or renting; Postal services; Subscriptions |
| Telecommunication |
| Aircraft renting or charter of aircraft |
| Automotive vehicle renting or charter |
| Taxi, toll and parking |
| Terrestrial, maritime and fluvial tickets |
| Watercraft renting or charter |
| Participation in course, talk or similar event |
Para fazer nossas comparações, precisaremos separar alguns conjuntos por regiões do país (vamos ver quem é bom de geografia!).
O operador | na função filter indica “ou”.
camara.NE = camara.sample%>%
filter(state == "RN" |state == "PE" |state == "BA" | state == "SE" |state == "AL" |state == "CE" | state == "PB" | state == "PI" |state == "MA")Acima, o conjunto camara.NE está recebendo os dados filtrados cuja coluna state esteja marcada com RN, ou PE, ou BA, ou…
Faça o mesmo para as outras regiões do país.
camara.NE = camara.sample%>%
filter(state == "RN" |state == "PE" |state == "BA" | state == "SE" |state == "AL" |state == "CE" | state == "PB" | state == "PI" |state == "MA")
camara.SE = camara.sample%>%
filter(state == "RJ" |state == "ES" |state == "MG" | state == "SP")
camara.N = camara.sample%>%
filter(state == "AM" |state == "PA" |state == "AP" | state == "RR" | state == "RO" | state == "AC" | state == "TO")
camara.CO = camara.sample%>%
filter(state == "MT" |state == "MS" |state == "GO" | state == "DF")
camara.S = camara.sample%>%
filter(state == "RS" |state == "SC" |state == "PR")Será que deputados do NE gastam mais com passagens aéreas que deputados do SE?
Primeiro precisamos ver a distribuição dos dados para passagens aéreas.
Como você faria esse gráfico?
Use unique() para ver o nome dado na tabela para os gastos com passagens aéreas.
Há dois gastos: usaremos ‘Flight ticket issue’.
Agora crie o conjunto de dados desses gastos.
NE.flight = camara.NE%>%
filter(subquota_description == "Flight ticket issue")
SE.flight = camara.SE%>%
filter(subquota_description == "Flight ticket issue")
boxplot(SE.flight$total_net_value, NE.flight$total_net_value)O boxplot mostra que os dados tem distribuição normal. O histograma confirma.
hist(SE.flight$total_net_value)O boxplot mostra que os dados tem distribuição normal. O histograma confirma.
hist(NE.flight$total_net_value)Que tipo de teste deve ser feito para comparar os gastos das duas regiões?
t.test(NE.flight$total_net_value, SE.flight$total_net_value, paired = F, alternative = "two.sided")##
## Welch Two Sample t-test
##
## data: NE.flight$total_net_value and SE.flight$total_net_value
## t = 0.5413, df = 694.87, p-value = 0.5885
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -58.26187 102.61541
## sample estimates:
## mean of x mean of y
## 409.6145 387.4378
Será que deputados do NE gastam mais com passagens aéreas que deputados do N?
Crie os dataframes necessários à análise
Faça os gráficos de diagnóstico de distribuição dos dados
Defina que teste usar
Faça o teste
N.flight = camara.N%>%
filter(subquota_description == "Flight ticket issue")
boxplot(N.flight$total_net_value, NE.flight$total_net_value)t.test(NE.flight$total_net_value, N.flight$total_net_value, paired = F, alternative = "two.sided")##
## Welch Two Sample t-test
##
## data: NE.flight$total_net_value and N.flight$total_net_value
## t = -2.5391, df = 232.3, p-value = 0.01177
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -315.15226 -39.75996
## sample estimates:
## mean of x mean of y
## 409.6145 587.0707
Será que deputados do RN e da PB gastam o mesmo com refeições?
Como fazer o conjunto de dados?
RN.meal = camara.sample%>%
filter(subquota_description == "Congressperson meal", state == "RN")
PB.meal = camara.sample%>%
filter(subquota_description == "Congressperson meal", state == "PB")boxplot(RN.meal$total_net_value, PB.meal$total_net_value)Como vimos na última aula, podemos fazer um teste de Wilcoxon
wilcox.test(RN.meal$total_net_value, PB.meal$total_net_value, alternative = "two.sided", paired = F)## Warning in wilcox.test.default(RN.meal$total_net_value,
## PB.meal$total_net_value, : cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: RN.meal$total_net_value and PB.meal$total_net_value
## W = 87, p-value = 0.1174
## alternative hypothesis: true location shift is not equal to 0
Às vezes precisamos que os dados tenham distribuição normal para fazer alguns testes específicos. Uma possibilidade, é fazer uma transformação logarítmica nos dados.
A função mutate() do dplyr cria uma nova coluna com o nome e os valores especificados. Aqui, criamos a coluna log.value que receberá os valores da transformação logarítmica da coluna total_net_value
Quando você rodar o código, verá que aparecerá uma nova coluna nos dataframse RN.meal e PB.meal descritos na área “Global Environment”.
RN.meal = RN.meal%>%
mutate(log.value = log(total_net_value))
PB.meal = PB.meal%>%
mutate(log.value = log(total_net_value))Podemos ver que os valores dessa nova coluna tem uma distribuição mais próxima a normal.
boxplot(RN.meal$log.value, PB.meal$log.value)Agora podemos fazer um teste-t
t.test(RN.meal$log.value, PB.meal$log.value, paired = F, alternative = "two.sided")##
## Welch Two Sample t-test
##
## data: RN.meal$log.value and PB.meal$log.value
## t = -1.4472, df = 27.63, p-value = 0.1591
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.2932263 0.2227894
## sample estimates:
## mean of x mean of y
## 3.797977 4.333196
Quando e/ou por que eu deveria transformar os meus dados?
Essa pergunta não tem uma resposta única ou fácil. No caso de dados que poderiam ser analisados pelo teste-t, há a possibilidade de usarmos um teste de Wilcoxon. Como há uma alternativa, não há um motivo muito forte para transformar o dados.
Há modelos mais complexos, que geralmente lidam com medidas repetidas, que pressupõem uma distribuição normal dos dados. Nesses casos, é comum lermos artigos que se baseiam na transformação dos dados.
O ideal é conhecer a literatura da sua área e conhecer bem os seus dados para tomar uma decisão.
É possível fazer quantos teste-t (ou Wilcoxon) eu quiser. Por exemplo, eu posso querer comparar os valores gastos pelos deputados de cada estado com bilhetes aéreas. Isso me daria quase 300 comparações!
Vamos fazer?!
Calma! O R automatiza esse processo para você através da função pairwise.t.test().
Digite ?pairwise.t.test para vermos os argumentos dessa função.
Vamos ter que criar um dataframe apenas com os dados de gastos com bilhetes aéreos.
camara.flight = camara.sample%>%
filter(subquota_description=="Flight ticket issue")pairwise.t.test(camara.flight$total_net_value, camara.flight$state, paired = F)##
## Pairwise comparisons using t tests with pooled SD
##
## data: camara.flight$total_net_value and camara.flight$state
##
## AC AL AM AP BA CE DF ES GO MA MG
## AC 1.000 - - - - - - - - - - -
## AL 1.000 1.000 - - - - - - - - - -
## AM 1.000 1.000 1.000 - - - - - - - - -
## AP 1.000 1.000 1.000 1.000 - - - - - - - -
## BA 1.000 1.000 1.000 1.000 1.000 - - - - - - -
## CE 1.000 1.000 1.000 1.000 1.000 1.000 - - - - - -
## DF 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - - - - -
## ES 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - - - -
## GO 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - - -
## MA 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - -
## MG 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -
## MS 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.538 0.489
## MT 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## PA 1.000 1.000 1.000 1.000 1.000 0.253 1.000 1.000 1.000 1.000 0.066 0.027
## PB 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## PE 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## PI 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## PR 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## RJ 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## RN 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## RO 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## RR 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## RS 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## SC 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## SE 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## SP 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## TO 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## MS MT PA PB PE PI PR RJ RN RO RR RS
## AC - - - - - - - - - - - -
## AL - - - - - - - - - - - -
## AM - - - - - - - - - - - -
## AP - - - - - - - - - - - -
## BA - - - - - - - - - - - -
## CE - - - - - - - - - - - -
## DF - - - - - - - - - - - -
## ES - - - - - - - - - - - -
## GO - - - - - - - - - - - -
## MA - - - - - - - - - - - -
## MG - - - - - - - - - - - -
## MS - - - - - - - - - - - -
## MT 1.000 - - - - - - - - - - -
## PA 1.000 1.000 - - - - - - - - - -
## PB 1.000 1.000 1.000 - - - - - - - - -
## PE 1.000 1.000 1.000 1.000 - - - - - - - -
## PI 1.000 1.000 1.000 1.000 1.000 - - - - - - -
## PR 1.000 1.000 1.000 1.000 1.000 1.000 - - - - - -
## RJ 0.417 1.000 0.018 1.000 1.000 1.000 1.000 - - - - -
## RN 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - - - -
## RO 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - - -
## RR 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 - -
## RS 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -
## SC 0.845 1.000 0.181 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## SE 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## SP 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## TO 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## SC SE SP
## AC - - -
## AL - - -
## AM - - -
## AP - - -
## BA - - -
## CE - - -
## DF - - -
## ES - - -
## GO - - -
## MA - - -
## MG - - -
## MS - - -
## MT - - -
## PA - - -
## PB - - -
## PE - - -
## PI - - -
## PR - - -
## RJ - - -
## RN - - -
## RO - - -
## RR - - -
## RS - - -
## SC - - -
## SE 1.000 - -
## SP 1.000 1.000 -
## TO 1.000 1.000 1.000
##
## P value adjustment method: holm
Vamos mudar a pergunta: suponha que você tem 20 bolas dentro de uma caixa, 19 pretas e 1 vermelha.
R: 0.05.
R: Aumenta! Se você tem 20 chances de tirar a bola vermelha, maior a chance e tirá-la!
Temos aqui o mesmo princípio das bolinhas vermelhas: https://xkcd.com/882/
É comum precisarmos fazer mais de uma comparação no mesmo conjunto de dados. Nesse caso, é preciso corrigir o nosso p-valor para que não aumentemos as nossas chances de falso positivo ou falso negativo.
Há várias opções para correção de p-valor. A escolha por uma delas deve levar em conta o quanto você quer controlar para falsos positivos e falsos negativos. Mais uma vez, leia artigos na sua área para ver o que as pessoas costumam usar. Não deixe de corrigir o seu p!
Confira o help da função p.adjust() para ver os métodos disponíveis no R.
Erro do Tipo I: falso positivo
Erro do Tipo II: falso negativo
Fonte: https://twitter.com/justinwolfers/status/666448547097677829
Será que deputados do NE e da SE gastam o mesmo com refeições?
Crie os conjuntos de dados
Faça um boxplot
SE.meal = camara.SE%>%
filter(subquota_description=="Congressperson meal")
NE.meal = camara.NE%>%
filter(subquota_description=="Congressperson meal")boxplot(NE.meal$total_net_value, SE.meal$total_net_value)Quase R$ 2000,00 em uma refeição?!
Vamos clicar na tabela para dar uma olhada nesses valores extremos.
Clique em total_net_value para organizar a lista em ordem crescente/decrescente.
Vamos excluir os valores extremos referentes a algo que provavelmente não é uma refeição.
NE.meal = camara.NE%>%
filter(subquota_description=="Congressperson meal", supplier != "VERDE VALE HOTEL S/A")Isso, claro, não muda o fato de a distribuição não ser normal
Distribuição não é normal: wilcoxon
wilcox.test(NE.meal$total_net_value, SE.meal$total_net_value, paired = F, alternative = "two.sided")##
## Wilcoxon rank sum test with continuity correction
##
## data: NE.meal$total_net_value and SE.meal$total_net_value
## W = 57994, p-value = 2.123e-06
## alternative hypothesis: true location shift is not equal to 0
median(NE.meal$total_net_value) ## [1] 60.5
median(SE.meal$total_net_value)## [1] 45.98
Os gastos dos deputados do Nordeste com alimentação foram significativamente maiores que os de deputados do Sudeste (mediana: 60,50 (n = 288) vs 45,98 (n = 331); W = 57994, p < 0,0001).
Imagine que você queira fazer os boxplot de gastos de refeição por estado. Primeiro você cria o conjunto de dados:
camara.meal = camara.sample%>%
filter(subquota_description == "Congressperson meal")E agora?
Fazer 28 gráficos, 1 a 1, dá trabalho. Felizmente, é só colocar “state” no eixo x do boxplot
ggplot(camara.meal, aes(y=total_net_value, x = state))+
geom_boxplot(outlier.size = 2, outlier.alpha = 0.3)Mas e se eu quiser fazer um histograma?
ggplot(camara.meal, aes(x=total_net_value))+
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Mas e se eu quiser fazer um histograma?
A função facet_wrap(~) faz isso, indicando que tudo que venha antes de ~ deverá ser agrupado pelo vetor indicado na sequência
ggplot(camara.meal, aes(x=total_net_value))+
geom_histogram()+
facet_wrap(~state)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Algumas observações:
geralmente o boxplot é uma maneira mais fácil de perceber diferenças “no olho”
o facet_wrap é uma opção para quando você quiser agrupar certos valores e mostrá-los em gráficos separados, porém lado a lado
Eu gostaria de saber as médias dos gastos em cada subquota. Qual seria um caminho possível para fazer isso?
agrupar os dados por subquota (group_by())
tirar a média de cada subquota (summarise())
camara.sample%>%
group_by(subquota_description)%>%
summarize(media = mean(total_net_value))## # A tibble: 22 x 2
## subquota_description media
## <fct> <dbl>
## 1 Aircraft renting or charter of aircraft 9281.
## 2 Automotive vehicle renting or charter 4453.
## 3 Automotive vehicle renting or watercraft charter 2583.
## 4 Congressperson meal 82.5
## 5 Consultancy, research and technical work 7114.
## 6 Flight ticket issue 434.
## 7 Flight tickets 2054.
## 8 Fuels and lubricants 609.
## 9 Locomotion, meal and lodging 1491.
## 10 Lodging, except for congressperson from Distrito Federal 409.
## # … with 12 more rows
Organizar a tabela em ordem decrescente pela média
camara.sample%>%
group_by(subquota_description)%>%
summarize(media = mean(total_net_value))%>%
arrange(desc(media))## # A tibble: 22 x 2
## subquota_description media
## <fct> <dbl>
## 1 Aircraft renting or charter of aircraft 9281.
## 2 Publicity of parliamentary activity 7142.
## 3 Consultancy, research and technical work 7114.
## 4 Automotive vehicle renting or charter 4453.
## 5 Participation in course, talk or similar event 3636.
## 6 Watercraft renting or charter 3011.
## 7 Automotive vehicle renting or watercraft charter 2583.
## 8 Flight tickets 2054.
## 9 Security service provided by specialized company 1938.
## 10 Locomotion, meal and lodging 1491.
## # … with 12 more rows
camara.sample%>%
group_by(congressperson_name)%>%
summarize(media = mean(total_net_value))%>%
arrange(desc(media))## # A tibble: 1,188 x 2
## congressperson_name media
## <fct> <dbl>
## 1 FERNANDO DE FABINHO 17325.
## 2 FRANCISCO RODRIGUES 12709.
## 3 ADÃO PRETTO 11995.
## 4 SABINO CASTELO BRANCO 11752.
## 5 AGNALDO MUNIZ 11408.
## 6 ARLINDO CHINAGLIA 11162.
## 7 JOSÉ MENDONÇA BEZERRA 9441.
## 8 SÁ 9231.
## 9 LUCIO VIEIRA LIMA 9087.
## 10 PAULO BAUER 8944.
## # … with 1,178 more rows
Grupo 1: analise se há diferença entre os gastos dos deputados do NE e SE para participação em congresso (“Participation in course, talk or similar event”). Apresente um boxplot comparativo, decida que teste usar e reporte o resultado. Faça também um gráfico comparativo do gasto por estado para cada uma das regiões, para ver se os estados diferem muito entre si.
Grupo 2: analise se há diferença entre os gastos dos deputados do SE e S para aluguel de carros (“Automotive vehicle renting or charter”). Apresente um boxplot comparativo, decida que teste usar e reporte o resultado. Faça também um gráfico comparativo do gasto por estado para cada uma das regiões, para ver se os estados diferem muito entre si.
Para gráficos
facet_wrap(~ x)Para manipulação de dados com dplyr
filter()
mutate()
groyp_by()
arrange()
desc()